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ABSTRACT 

We develop a model for the growth of dark matter halos and use it to study their 
evolved density profiles. In this model, halos are spherical and form by quiescent 
accretion of matter in clumps, called satellites. The halo mass as a function of redshift 
is given by the mass of the most massive progenitor, and is determined from Monte- 
Carlo realizations of the merger-history tree. Inside the halo, satellites move under 
the action of the gravitational force of the halo and a dynamical friction drag force. 
The associated equation of motion is solved numerically. The energy lost to dynamical 
friction is transferred to the halo in the form of kinetic energy. As they sink into 
the halo, satellites continually lose matter as a result of tidal stripping. The stripped 
matter moves inside the halo free of dynamical friction. The evolved density profiles 
are steeper than those obtained by assuming that, once they have been accreted onto 
the parent halo, satellites remain at a fixed distance from the halo center. We find 
that the final density profile depends mainly on the rate of infall of matter onto the 
halo. This, in turn, depends on the initial fluctuation field as well as on cosmology. 
For mass scales where the effective spectral index of the initial density field is less 
than —1, the model predicts a profile which can only approximately be matched by 
the one parameter family of curves suggested by Navarro, Frenk and White (1997). 
For scale-free power-spectra with initial slope n, the density profile within about 1% 
of the virial radius is p oc r _/5 , with 3(3 + n)/(5 + n) < (3 < 3(3 + n)/(4 + n). 
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1 INTRODUCTION 

A long standing question in the dynamics of cosmological gravitating systems is: how well do dynamically evolved systems 
retain memory of their initial conditions? On large scales, where the evolution is still in the the weakly non-linear regime, 
the growing mode of the initial density fluctuations can, in principle, be fully recovered, if the present velocity or density 
field is given (e.g., Peebles 1989, Nusser & Dekel 1992). On small scales where shell crossing has occurred and "virialised" 
objects (halos) have formed, the situation is less clear. Using techniques of statistical mechanics, Lynden-Bell (1967) showed 
that starting from a general initial state, a gravitating system could, via a process he termed "violent relaxation", reach a 
quasi-equilibrium state which is almost independent of the initial conditions. Therefore, under the restricted conditions he 
assumed in his analysis, a gravitating system develops a "universal" density profile independent of the initial conditions. 

The assumptions underlying Lynden-Bell's analysis are hard to justify in the hierarchical scenario for the formation 
of structure in an expanding universe. Nevertheless, high resolution cosmological N-body simulations of the gravitational 
clustering of collisionless particles from hierarchical initial conditions strongly suggest that dark matter halos do indeed 
develop a universal final density profile (e.g. Dubinski & Calberg 1991, Lemson 1995, Navarro Frenk & White 1995, 1996, 
Cole & Lacey 1996, Moore et. al. 1997). Navarro, Frenk & White (1995, 1996; hereafter NFW) found that density profiles of 
halos of different masses in a variety of cosmological models could be fit with the following one parameter functional fit 
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where pb is the background density and r v is the virial radius of the halo, defined as the radius within which the average 
density is 178 times that of the background. With this definition for r v , the parameter 5n can be expressed in terms of the 
concentration parameter c; thus, c is the only free parameter in the fit ([j]). It has been argued by NFW that c is directly 
related to the formation time of a given halo. The claim that density profiles could be fitted with a one parameter functional 
form, has recently been challenged by Klypin et. al. (1998) who found, using N-body simulations of CDM-like models, that 
the scatter about a one parameter fit is substantial, indicating that the structure of halo density profiles involves more than 
just one physical parameter. 

Because of the lack of a general analytic technique for following the detailed evolution of a system from general initial 
conditions, most analytic work has focused on studying the evolution of isolated spherical systems (Gunn & Gott 1972). The 
collapse of a spherical density perturbation of a self-similar form, S oc r~ m , in an otherwise flat universe, yields the density 
profile p oc r ' 3m/{1+m) for m > 2, and p oc r~ 2 for m < 2 (Filmore & Goldreich 1984, Bertschinger 1985). This implies 
that in highly non-linear systems, full information about the initial distribution is preserved only for a special class of initial 
conditions. Strictly speaking, this result is valid only for the case of purely radial collapse. In principle, non-radial motions 
can prevent particles with large turnaround radii from sinking to the inner regions of the collapse and forming an r~ 2 profile. 
If particles were assigned angular momenta in a self-similar way, then the density profile above is expected to be valid when 
m < 2 as well (White & Zaritsky 1992). 

These special spherical solutions were first related to the formation of dark matter halos from initial gaussian density 
fields by Hoffman & Shaham (1985) (also see Hoffman 1988). They noted that the mean shape of high peaks in gaussian 
fields is S oc r~^ n+: ^ where n is the index of the initial power spectrum (Dekel 1981, Bardeen et al. 1986). Therefore, they 
argued that setting m = n + 3 makes the spherical solutions relevant in the cosmological context. There are two caveats to 
this argument. First, it assumes that the accreting matter is not clumpy. In contrast, in hierarchical scenarios for structure 
formation, halos grow in mass by mergers with other, typically less massive, halos (cf. Lacey & Cole 1994; Lemson 1995; Syer 
& White 1996). This means that the effects of dynamical friction and tidal stripping may well be important in determining the 
final structures of halos in hierarchical models. Second, one is usually interested in the density profiles of halos of a given mass 
today; these may have had, as their seeds, initial density peaks of various heights and shapes. Moreover, strictly speaking, 
the Hoffman & Shaham (1985) solution applies only in regions outside the virial radius of a dark matter halo. In fact, Syer & 
White (1996) argue that the profile must scale as r -13 ^ with /few = 3m/(2 + m) in the inner regions of dark matter halos, 
and this solution differs from the r -13 ^ , with /3hs = — 3m/(l + m), scaling proposed by Hoffmann & Shaham (1985). 

Since understanding the density profiles of halos is necessary to relate structure formation theory to the observed structure 
of galaxies, we aim here to develop a more detailed model of the process of halo formation. We will continue to assume spherical 
symmetry, and will formulate a semi-analytic scheme for modelling the evolution and the structure of halos more realistically. 
The first goal of this paper is to develop and test a method for tracing the mass of a given halo back in time that is consistent 
with the hierarchical clustering scenario. We describe this method briefly in section 2.1 (details are discussed in the Appendix). 
Our second goal is to provide a simple model, which includes the effects of dynamical friction and tidal stripping, for the 
motion of satellites once they are inside a halo, and to then use this model to estimate the halo structure that results. In 
section 2.2 we describe how to treat the motion of satellites once they are inside the halo. The equation of motion must be 
solved numerically and, in section 3, we present the results of applying this scheme to simulate the evolution of a variety of 
halos. We study halos of various masses that form from initially scale free gaussian fluctuation fields as well as fields with the 
CDM power spectrum. In section 4 we conclude with a discussion of our results and argue that the concentration scale c in 
the NFW profile, if fundamental, cannot directly be related to the formation time and must reflect an interplay between the 
formation time and other — yet unknown — physical effects, which determine the halo profile. 



2 THE MODEL 

We will compute the evolved density profile of a halo which has mass Mo at the final redshift zo in two steps. First, we devise 
a scheme for tracing back, to some early time, the mass of the largest progenitor subclump of the halo. The rate of change of 
the mass of the most massive progenitor (MMP) gives the accretion rate of matter onto the already existing halo. This rate 
is sensitive to the background cosmology, and to the shape of the initial power spectrum. Our second step is to formulate a 
dynamical prescription for following the evolution of the accreted matter inside the halo. 
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2.1 Halo Mass Growth Rate 

At any given redshift z, define the current virial radius, r v , of a halo as the radius within which the average density is 178pb, 
where pt, is the mean density in the universe at that time. At z, the current halo mass is defined as the mass within the 
current virial radius. Ofcourse, the current virial radius and mass are smaller than their final values. Two key ingredients 
in our model are (i) how the mass within the virial radius grows with time, and (ii) how the accreting mass is distributed 
among satellites. To provide these two ingredients, one must trace the full merger history of a given halo back in time. That 
is, we need to know how the mass Mo was partitioned into progenitor halos at any earlier time Z\ > Zq. For Gaussian initial 
conditions, merger history trees of halos, each of which has the same mass Mo at zo, can differ substantially. Therefore, we 
need to list all possible merger trees, and we need to compute the probability that each occurs. Except for Poisson initial 
conditions (Sheth 1996), this probability has not been computed analytically. Hence, we adopt a short-cut. We trace only the 
most massive progenitor of the most massive progenitor, and so on, back in time. 

We do this as follows. We assume that at z\, the average number of progenitors of an (Mo,zo)-halo that have mass 
between M\ and (Mi + dMi) can be approximated by 



N(M 1 , Zl \M ,z )dM 1 = (^f) * ^-^ exp 
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(Bower 1991, Lacey & Cole 1993), where \fSl is the rms density fluctuation in a Top-Hat window function of radius 
(3M 1 /4tt) 1/3 . Lacey & Cole (1994) show that this expression is in good agreement with what happens in numerical sim- 
ulations (but see Tormen 1998). The second equality defines /(1|0), which is the fraction of the mass of Mo that, at zi, was in 
objects with mass Mi. This means, of course, that the integral of /(1|0) over the range < Mi < Mo equals unity. However, 
as (zi — zo) — * 0, the integral of N(Mi, zi|Mo, zo) over the range (Mo/2) to Mo also approaches unity. Since an Mo halo 
may have at most one progenitor with mass in this range, TV (1 1 0) can be interpreted as the probability that an (Mo, zo)-halo 
had an (Mi, 21) progenitor subhalo (cf. Lacey & Cole 1993), with (Mi/Mo) restricted to the range between one half and 
unity, a short time earlier. Consequently, in the limit of small redshift intervals, we choose the mass of the MMP according 
to equation (|^). We then replace Mo with the chosen value of Mi, and z a with z\ > zo, and iterate until Mi, the mass of the 
MMP is as small as desired. In the Appendix we argue that this scheme provides MMP histories that are similar to those 
which occur in numerical simulations. 

In what follows, the ensemble average of the accretion rate will be useful. Before computing it, the following question 
arises: should we use curves of M(z) by averaging the masses of all realisations of merger histories in the mass direction given 
the redshift, or should we use M(z) obtained by averaging the redshift at a given mass? Fortunately, the two ways of averaging 
lead to almost identical curves of M(z). Figure 1 shows curves of M(z) for halos having Mo/M* =1,3 and 10, for scale free 
power spectra with n = 0, —1 and —2. Also shown are results for halos having M = 10 11 , 10 13 and 1O 15 M0 in the standard 
CDM model, normalised such that the linear rms value of density fluctuations in a spherical window of radius 8h~ Mpc is 
us = 0.6. The solid lines show curves obtained by averaging the mass at a given redshift, and the dashed lines show curves 
obtained by averaging the redshift at a given mass. The error bars are 1-cr deviations from the mean M(z). The dashed and 
the solid lines corresponding to each case almost overlap. This figure also shows the well known fact that, on average, lower 
mass halos formed at an earlier epoch. In what follows we have arbitrarily chosen to work with M(z) obtained by averaging 
the merger histories in the mass direction. ^] 

2.2 Density profiles from stable clustering 

One way to convert the mass growth curves M(z) to density profiles is as follows. We assume that the mass M(z~ dz) — M(z) 
accreted during dz is accreted smoothly onto the MMP. This means that we are ignoring the fact that this mass may actually 
be divided among subclumps. Assume that a newly accreted patch of matter oscillates around the center of the halo with 
mean radius equal to the virial radius of the halo at the redshift at which it joined the halo. This is the stable clustering 
assumption, since it is equivalent to assuming that there is no net in-flow of matter in the virialised region. It is motivated by 

* When the mass of the MMP is greater than half of Mo, then M(z) can be computed directly from the formation time argument given 
in the Appendix. The mean value of the scaled variable there, ajf , is 

. fP (. f )d. f = yiQ + _L), (n = 0) 

from which the mean value of z given M(z)/Mq is easily obtained. 
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Figure 1. Curves of the average M{z)/Mq for power law spectra with n = —2, -1 and 0, as indicated in the plot. The solid curves, in 
each panel, correspond, from top to bottom, to M/M„ = 0.1, 1 and 10 obtained by averaging in the mass direction. The dotted lines 
in the middle panel lines are la deviations about the mean M(z) for 1M, halo. The dashed line in the middle panel is M(z) for 1M, 
obtained by averaging in the redshift direction. 



the analytic solutions for collapse of density peaks with pure power law initial profiles (Filmore & Goldreich 1984; Bertschinger 
1985). Recall that, in these solutions, the mean density within a radius r of the halo is a multiple (here we use 178) of the 
background density of the universe at the time when the virial radius of the halo was r. Indeed, in this approximation, our 
only change to the Hoffman & Shaham (1985) study is our method for computing the accretion rate. A similar analysis was 
also done by Avila- Reese, Firmani & Hernandez (1997) who used the prescription of Lacey & Cole (1993) and Eisenstein & 
Loeb (1996) to compute M(z). 

Ideally, one would like to use stable clustering to compute final profiles for each individual merger history and then 
average the result over many realisations. In practice this is time consuming. We show below (Figure 4) that the evolved 
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Figure 2. Density profiles obtained from curves of M(z) under the assumption of stable clustering. The dotted curves are the NFW 
best fits to the density profiles. 



density profile computed from the ensemble averaged M(z) is a good approximation to the ensemble average of the individual 
evolved profiles. To do so, we have used the stable clustering assumption to compute profiles of 10 random realisations of the 
merger tree of an M* halo for a scale free spectrum with n = — 1. The average of these profiles is the dashed curve in figure 4. 
The profiles from six individual merger histories are shown by the light solid curves. The thick solid line is the profile obtained 
from the average M(z) shown in figure 1. We conclude from figure 4 that, at least under the assumption of stable clustering, 
working with the average M(z) does not introduce systematic biases. 

The density profiles generated in this way from the ensemble averaged M(z) curves of Fig. 1 are shown in Fig. 2. The dotted 
lines are fits of the form NFW form given in ([l]). Figure 3 shows the associated circular velocity profiles: V^{r) oc GM(r)/r. 
Except when n = 0, the individual density and circular velocity profiles are reasonably well described by the NFW parametric 
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Figure 3. The same as the previous figure but showing rotation curves instead of density profiles. 



form. Notice that less massive halos are more centrally concentrated. The trend of steeper profiles for larger values of the 
spectral index, n, is also seen in the simulations (e.g. Cole & Lacey 1996). There is, however, a fundamental difference between 
the profiles obtained here and those in the simulations. For example, in the CDM N-body simulations, profiles of halos with 

masses ~ 10 11 M Q ~ 1O 15 M almost coincide over the range 0.1 > r/r v < 1 (cf. figure 4 in NFW 1995). On the other 

hand, the corresponding profiles in figure 2 substantially deviate from each other for the same range. Thus, the concentration 
parameter c depends differently on mass in our models than it does in the simulations. NFW argue that c is determined by 
the formation time of halos. Our detailed models of the accretion history incorporate this formation time explicitly, yet, they 
are unable to reproduce the trends of c with mass that are seen in simulations. This suggests that, if simulated halos are 
indeed well described by the one parameter c, then this parameter must depend on more than just the formation time; c must 
depend on some additional physics. 
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Figure 4. Density profiles from 6 random realisations (points) of M(z) for n = — 1 and 1M* halo. Plotted also is the average profile 
(heavy dashed) of 10 random realisations and the profile (heavy solid) corresponding to the average M(z). The profiles were computed 
under the stable clustering assumption. 



The stable clustering assumption is not expected to hold in realistic collapse configurations. This is because the accreting 
matter is, in fact, distributed into bound clumps (satellites). These will suffer from dynamical friction as they fall into the halo. 
This causes a transfer of energy from the satellites to the halo which is neglected in the stable clustering collapse considered 
above. Moreover, the tidal field of the halo continuously prunes infalling satellites, and this is also ignored in the smooth 
collapse scenario. In the next section we describe how to compute evolved density profiles which include these effects. 
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2.3 Beyond stable clustering: dynamical friction, tidal stripping and halo heating 

In this section we formulate a dynamical model of the motion of satellites inside halos. The treatment here does not require 
any assumptions about the matter accretion rate, provided major mergers do not occur. Indeed mergers of two halos of similar 
mass invalidate the dynamical prescription presented here. 

We will write the equations of motion of satellites inside a halo under the approximation that the halo is spherical. Then, 
the gravitational field of the halo is that of a spherical mass distribution. In addition to the gravitational field of the halo, 
satellites suffer from a dynamical friction drag force (Chandrasekhar 1943). Therefore, the equations governing the motion of 
a satellite within a halo are, 

d 2 r . . dr , , , . r ,„, 

— = -u A1 (r,t)--M h (r,t)^, (3) 

where Mh(r, t) is the mass of the halo within radius r at time t and Wdf is the inverse of the dynamical friction time.Q Ostriker 
& Turner (1979) derived a similar expression. For a satellite of mass M s moving in a medium of density, p, and velocity 
dispersion, a, the quantity w<jf j can be approximated by 

Wdf (r,i)=AfM s 4r4 (4) 

(Binney & Tremaine 1987), where N is a numerical factor. Here cr(r,t) and p(r,t) represent the average velocity dispersion 
and density within r. 

The gravitational field of the halo prunes the satellites as they sink to the center of the halo. This causes a gradual 
decrease in the mass of the satellite, M B , which should be taken into account in the equations of motion (^). We assume that 
a satellite at a distance r from the center of the halo retains only matter which lies inside a certain radius, the tidal radius 
r t , from its center. This tidal radius is approximately that radius at which the gravitational field of the satellite equals the 
change in the gravitational field of the main halo times rt. That is, 

d_(M^\ M s 
dr 



fMh\ M s 



which, to first order in r t /r, reduces to 
M.(r t ) _ M h (r) 



(6) 



Thus, the average density of the satellite within rt equals the average density of the halo within r. Given a guess about the 
density profile of the satellite when it crosses the virial radius of the halo, and assuming that the profile within r t does not 
change with time, the relation (|^) determines the mass of the satellite as a function of r. 

Here we assume that satellites are spherical and have 1/r 2 density profiles, i.e., M s (r t ) oc r t . In this case, as a satellite of 
mass M sv passes through a halo of density p v , it loses mass according to 

. 1/2 

Pv 



where M s (r) is the mass of the satellites mass when it is at radius r and ph(r) is the average density of the halo within that 
radius. Substituting this in (Q), and using the fact that Mh(r)/r oc o 2 (r), we find that 

( p v \ 1/2 M hv 

Wdf = Q UJ Wdf)' (8) 

where a = \ // 3tvAT M BV / Mhv , and Mhv is the mass of the halo when the satellite crossed its virial radius. The factor a is scaled 
so that it is the ratio of the circular orbit time to the dynamical friction time at the time when the satellite joined the halo. 

Matter stripped from satellites is said to be 'evaporated'. Initially, evaporated matter is assigned the velocity and position 
of the satellite from which it was stripped. Thereafter, it moves under the influence of the gravitational field of the halo free 
of dynamical friction. Of course, the satellite, now with reduced mass, continues to suffer from dynamical friction. 

Next, we need to determine the fate of the energy that is lost by satellites as a result of the dynamical friction term in 
the equation of motion (^). A satellite loses energy at the rate 

d(ig/M.) _ /dr\ 2 



t We work in units in which (3=1. 
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(cf. Landau & Lifshitz 1960). It is unclear how this lost energy is redistributed within the halo. Since the bulk of the halo's 
matter contributes to the dynamical friction drag force which is exerted on the satellites, it seems reasonable to assume that 
the lost energy is distributed equally among the halo matter. That is, the energy deposited to a patch of evaporated matter 
is assumed to be proportional to its mass. 



2.3.1 Numerical Scheme 

The equations of motion above can be solved numerically, provided that we make a number of additional assumptions which 
we describe below. 

Our numerical scheme simulates a halo made of particles of identical mass. Some of these particles are assigned to satellites 
and suffer from dynamical friction; the other particles move under the influence of the gravitational field of the halo only 
and are called halo particles. The code treats these two types of particles differently. Since a satellite particle may be tidally 
stripped from its parent, we need some way of deciding when this happens; once stripped the particle becomes a halo particle, 
and can never be captured by satellites. Dynamical friction provides a mechanism for the transfer of energy from satellites 
to the halo; the code allows for this by assuming that only halo particles can absorb this energy. We describe how both these 
effects are included in more detail below. 

We follow the evolution of a halo starting at an early time at which the halo mass was a tiny fraction of its final mass. 
Given the growth of the mass of the halo with time, we can compute the time at which each particle will cross the virial 
radius of the halo. At any given time, we simulate only the trajectories of those particles which are within the current virial 
radius. We assume that, when first accreted, particles move on circular orbits. 

Since we are only simulating particles within the current virial radius, some of these particles may have already been 
stripped from their parent satellites. Here we assume that about half the particles joining the halo in any given time step 
are associated with satellites and the other half are halo particles. We need to specify the parameter a for all those particles 
which belong to satellites. Notice that equations Q and (^) involve the mass of a satellite only through the parameter a 
which determines the dynamical friction time scale. This simplifies the problem considerably. In the numerical scheme, each 
satellite's particle is assigned a value of a drawn from a uniform distribution with a minimum value of zero and some maximum 
value, a max , which determines the total energy lost by the satellites. 

This choice implies that, at any time, the average satellite mass is proportional to the mass of the MMP. This is in the 
spirit of hierarchical clustering and is motivated by results from N-body simulations (Tormen 1997) . Now we need to decide 
which particles are stripped, and when. At any given time, each satellite particle is assigned a probability of remaining bound 
to its parent. According to equation (Q) this probability is proportional to the ratio of the average density within the particle's 
radius at the current, to that at the previous time step. 

Finally, in order to avoid numerical instabilities, the gravitational force acting on a particle at radius r, is softened 
according to 

F ^ = -j^f^jm (10) 

where m is the particle's mass, Ni is the number of particles with r < rt and e is the force softening parameter. For a closed 
system of particles without dynamical friction this force law conserves the total "energy" 

" f i mNi 



Etot = m 



2 , .2^/2 



(r? + e 



(11) 



(White 1983). The factor uidi, which determines the amplitude of the dynamical friction drag force, is evaluated according to 

1/2 M hvi 

* = ^ (12) 

where Mhvi and p v i are the mass and the average density of the halo within the virial radius at the time the satellite was 
accreted onto the halo. The sorting routine INDEXX in Press et. al. (1992) reduces the evaluation of the gravitational and 
dynamical friction forces to an TV log 2 TV process, where N is the total number of particles in the halo at any given time step. 
Once the forces have been computed, a leapfrog time integration scheme is used to move particles from the current time 
step to the next. In each time step the energy lost by satellite particles is computed using (ph. The satellites lost energy is 
distributed equally among the halo particles in the form of kinetic energy. Namely, halo particles absorb the energy lost by 
satellites by receiving "kicks" in random directions. 

Increasing a ma x increases the total amount of energy lost by the satellites. Unfortunately, the code is unstable for values 
of Qmax which lead to energy loss exceeding a significant fraction of the absolute value of the total energy. We will see that 
even a small amount of heating is sufficient to significantly change the halo density profile in the inner region relative to the 
stable clustering profiles. 
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Figure 5. Orbits of satellites with various values of a. The stable clustering density profile of a CDM 10 1 5Mq was used to integrate 
the equations. The bottom panel shows an orbit if satellites did not lose mass as a result of tidal stripping. 

3 RESULTS 

We ran simulations with 80000 particles in the final halo and with 30000 equally-spaced time steps. In all the runs we set the 
softening parameter e = 0.0005 r v , where r v is the virial radius of the halo at the final time. For convenience we adjust the 
numerical factors in the equation such that a is the ratio of the dynamical friction time to the circular orbit time when the 
satellite just joined the halo. Following each time step, the satellites' lost energy is reassigned to particles in the halo. 

It is interesting to inspect the effects of dynamical friction and tidal stripping. In figure 5, we plot trajectories of satellites 
moving in a halo with a density profile given by that of a 10 15 Mq CDM halo as obtained assuming stable clustering and 
shown in Fig. 2. The satellites were initially placed on circular orbits in the XY plane. The orbits were integrated with a — 3 
and 1/2 respectively. For comparison a trajectory with no tidal stripping is also shown. 
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Table 1. Slopes of density profiles in the region It)- 2 < r/r v < 10 _1 . For power law power spectra, M\ , M2 and M3 correspond to halos 
of masses 0.1, 1 and 10 in units of M*. For CDM, they correspond to 10 11 , 10 13 and 1O 15 M0. Top and bottom values in each mass row 
are for the stable clustering and the semi-analytic dynamical model, respectively. 
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Fig. 6 shows density profiles of halos at the final time. Results are shown for halos with M/M* = 0.1, 1 and 10 for the 
scale free power spectra and M = 10 11 , 10 13 and 1O 15 M0 for the standard CDM model. The dotted lines show the best fit 
NFW profiles. The curves were computed assuming the satellites' lost energy was redistributed in the form of kinetic energy. 
The profiles in Fig. 6 were computed assuming a max = 1/2. For the CDM model, the fraction of the total energy lost by 
satellites for this value of a ma x was 0.23, 0.11 and 0.04, for halos of mass 10 11 , 10 13 and 10 15 M Q , respectively. Fig. 7 shows 
the circular velocity profiles corresponding to the density profiles shown in Fig. 6. 

Relative to the profiles obtained assuming stable clustering, these density profiles are steeper, at least in the inner 
regions. This is because the halo is heated by the infalling satellites. Recall that, in our prescription, energy lost by the 
satellites produces random motions among the halo particles. If the density profile is shallower that r~ 2 , then the number of 
particles crossing a radius r inward, as a result of the added random motions, is larger than the number of particles crossing 
outward. We have tried various other heating recipes such as redistributing only half the lost energy in the form of kinetic 
energy while the other half in the form of potential energy (by appropriately rescaling the radii of the halo particles). The 
final profiles in all the heating recipes we adopted were similar but not identical; they all were steeper than the corresponding 
profiles obtained from stable clustering. 

Although the density profiles for n > —2, are not very well described by NFW fits, the density profiles for n = — 2 and 
CDM are. The corresponding circular velocity profiles, plotted in Fig. 7, show significant deviations from the NFW fits. The 
slopes of the density profiles in the inner regions r < 0.1r v are listed in table 1. The top and bottom numbers in each row 
correspond to the profiles obtained from stable clustering (Fig. 2) and from the profiles in Fig. 6, respectively. It is interesting 
that the slope for the halo of mass 1O 15 M in the CDM model agrees with the result of the high resolution N-body simulation 
of Moore et. al. (1997). 



4 SUMMARY AND CONCLUSIONS 

Assuming gaussian initial conditions, we addressed two related issues. These were: a) how to generate random realisations of 
the growth history of a halo whose mass at the present time is given, and b) how the density profile of a halo can be related 
to its merging history. 

The way we generate random realisations of halo histories is, in principle, different from that previously presented in the 
literature. Lacey & Cole (1993) and Eisenstein and Loeb (1996) choose to work with the function /(Mi, zi|Mo, 20) which is 
the probability that a mass element of halo Mo at zo is incorporated in a halo of mass Mi at z\ > zo- Instead, we work with 
N(Mi, 2i|Mo, zo) — (Mo/Mi)f(Mi,zi\Mo,zo), which can be interpreted as the average number of progenitors of Mo that 
have mass (Mi) at z\. We have checked that these previous methods lead to results that are similar to ours, and we show in 
the Appendix why this is so. 

Given the mass growth rate we presented results using two models for computing the evolved density profiles. The first 
was based on the stable clustering assumption in which a mass element is effectively assumed to remain at the virial radius 
when it joined the halo. The second invoked a more detailed description of the motion of satellites once they have crossed the 
virial radius of the halo, incorporating tidal stripping, dynamical friction and heating of halo particles by infalling satellites. 
We found that including these effects leads to steeper profiles in the inner regions than in the stable clustering case. For initial 
power-spectra P(k) oc k n , we find that the slope in the inner regions depends on n, and also on the final mass. Table 1 shows 
that the inner profile is p oc r _/3 , with /3sw < P < /fes- For a given power spectrum, the dependence of density profiles on 
the halo mass in our models is stronger than that seen in the simulations (e.g. NFW). In particular, for the CDM power 
spectrum, the relative difference between the profiles of halos of masses 10 11 M Q and 1O 15 M is about an order of magnitude at 
r — 0.1r v . On the other hand, the corresponding profiles in the simulations almost coincide for 0.1 < r/r v < 1. This suggests 

© 0000 RAS, MNRAS 000, 000-000 



12 Nusser & Sheth 




-2 -1 -2 -1 

log 10 [r/r v ] log 10 [r/r v ] 



Figure 6. Density profiles obtained from the dynamical model using curves of M(z) from Fig. 2. The dotted curves show the best NFW 
fits. 



that the NFW argument that the concentration parameter is entirely determined by the formation time is inadequate. It is 
likely that effects, in addition to those included in our models, must play a significant role in determining halo density profiles. 
For example, in their study of fully three-dimensional collapse, Huss et al. (1998) suggest that a bar-related instability may 
play a major role in determining the final profile. 

Using our simplified models, we hoped to shed some light on the current — arguably contradictory — results of various 
recent N-body simulations (Navarro Frenk & White 1997, Kravtsov et. al. 1997, Moore et. al. 1997) of halo density profiles. 
Based on these models, the general conclusion is that a one parameter form, such as the fit proposed by Navarro, Frenk & 
White (1997), is inadequate to describe the profiles for all radii and all power spectra. 
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Figure 7. The same as the previous figure but showing circular velocity profiles. 
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APPENDIX A: THE ALGORITHM 

This Appendix describes our MMP generating algorithm in detail. First we compare our algorithm with others that have been 
proposed in the literature. Then we show that, in at least one respect, it is able to reproduce results measured in numerical 
simulations of hierarchical clustering. This motivates further study: we show that, for our algorithm, the distribution of 
formation times and the distribution of formation masses can be computed analytically. 



Al Comparison with previous work 

Recall that the MMP is followed through a succession of small timesteps. At each step, the mass decreases according to 
equation (^). This means that the mass of the MMP is never smaller than half the mass of the parent halo at the previous 
time step. It is hard to justify this assumption on general theoretical grounds (but see Sheth & Pitman 1997 for a self- 
consistent, binary split model of the merger tree described by Sheth 1996). Nevertheless, this assumption is in agreement with 
results from N-body simulations; Fig. 5 in Tormen et al. (1997) shows that the masses of the MMPs in their simulations do 
not decrease by more than a factor of two at each small time step. 

Our algorithm differs slightly from an MMP generating algorithm that was suggested by Lacey & Cole (1993), and used by 
Eisenstein & Loeb (1996). They split the MMP according to /(1|0) = (Mi/M ) JV(1|Q), with no restriction on Mi/M [recall 
that the integral of /(1|0) over the range < Mi /Mo < 1 equals unity as required], but they then choose the larger of the two 
pieces. This means that their MMP has mass Mi with probability (Mi/M ) N(Mi\M ) + (Mo - Mi)/M N(M - Mi|M ), 
and Mi > M /2 always. In the limit (zi - z ) -> 0, N(Mi\M ) w N(M - Mi|M ), for a large range of Mi (cf. Fig. 7 in Sheth 
& Pitman 1997), so this is similar to using A r (l|0) for the split rule. The approximation of equality is good for white-noise 
initial conditions, but it becomes worse for power spectra with n < 0. 

Let t{ denote the formation time, defined as the time when the mass m of the MMP is a fraction / of the halo at the 
present time. Provided / = m/Mo > 0.5, this distribution can also be computed using the arguments presented in section 2.5.2 
of Lacey & Cole (1993). For white noise initial conditions (initial power spectra with slope n = 0), this distribution can be 
computed analytically. If 

u)t = S c0 (zi - zq)/^J cr?n - cr 2 Mo , 
then, provided / > 0.5, 

Pfcr) - {j - l) *<* (^) + (?-j)]fl »P (-4) , (» - 0)- (AD 
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Figure Al. The distribution of the formation time of halos, expressed in terms of the variable u> = z (M* /M) 1 / 3 / \J (2 2 / 3 — 1) for a 
power-law spectrum with slope n = —1. The dotted, solid and dashed curves correspond to M/M* = 0.1, 1 and 10, respectively. The 
points show equation (2.33) of Lacey &; Cole (1993) for n = 0. 



When / = 0.5 this reduces to equation (2.33) of Lacey & Cole (1993). For more general initial conditions the answer must 
be computed numerically. For all / > 0.5, the distribution for other values of n is very similar to that for n = (Fig. 7 
of Lacey & Cole 1993 shows this for / = 0.5). At least when / = 0.5, this distribution for the formation times fits the 
distribution measured in N-body simulations well (Lacey & Cole 1994). Since it is trivial to compute this distribution from 
our Monte-Carlo MMP generating algorithm (simply follow the MMP till its mass drops to / times the initial mass), it is 
important to check that our MMP generating algorithm is able to reproduce this distribution. 

As a test of our MMP history scheme, Fig. Al compares the distribution of formation times p(tf) obtained from our 



Monte-Carlo realizations with the analytic prediction given by (Al). The distributions are plotted in terms of the scaled 
variable cj t , and we have set / = 0.5. Thus, (t /t { ) 2/3 = 1 + ^V2 Q - 1(M /M*)" q/2 where a = (n + 3)/3 and n is the slope of 
the initial power spectrum. The figure shows results for n = — 1, although our Monte-Carlo distributions for other values of 
n are also in good agreement with the analytic prediction, as are the distributions for other values of / > 0.5. The scatter of 
the Monte-Carlo distributions about the analytic prediction is also comparable to that seen in the N-body simulations. We 
use this fact to argue that our scheme provides MMP histories that are similar to those in the numerical simulations. 



Fig. Al shows the formation time distribution obtained directly from the numerical Monte-Carlo code. Since the distri- 



bution is in reasonable agreement with that found in N-body simulations, it may be worth studying our algorithm further. 
Below, we show how to derive explicit analytic expressions for the distribution of formation times and masses. 



© 0000 RAS, MNRAS 000, 000-000 



16 Nusser & Sheth 



A2 Analytic calculation of the formation time distribution 

Suppose a random variable Y is drawn from a Gaussian distribution with mean zero and variance \. Then Y 2 has a 
Gamma(|,x) distribution, and X — Y~ 2 is said to have the Inverse Gamma(|,x) distribution. That is, 



Consider a sequence of random variables Xi,Xt, with each of the XiS chosen from an Inverse Gamma(|,x) distribution 
independently of the others. Define 

n 
i = l 

It is easy to verify that the distribution of S n is Inverse Gamma(i, n 2 x), so that 

m sj { n2 x/ x ) 1/2 ( n2 x\ da; , . t . 

p n (S„ =x)dx= - — ^=d — exp — — - — . (A4) 

\ 2x x 



These expressions will be useful below. 

Consider a halo that has mass Mo at zq. Suppose that at Z\ = zo + Az, the MMP of this halo has mass < Mi < Mo with 
probability f(Mi,zi\Mo,Zo), where /(1|0) is given by equation (Q). Note that, in the main text, we use a different function 
for the MMP probability. We use this form here to illustrate the way the calculation is done, and show the result of using our 
rule later. Then 

(Si -So) = (S c0 Az) 2 /y 2 , (A5) 

where y is drawn from a Gaussian distribution with zero mean and unit variance, so y 2 is drawn from a Gamma(|,l) 
distribution, which means that 1/y 2 is drawn from an Inverse Gamma(|, 1) distribution. Thus, the distribution of (Si — So) 
is Inverse Gamma(|,x), with \ ~ (^co A2) 2 . 

Suppose that, at 22 = 21 + Az, the MMP of Mi is obtained by drawing M 2 from /(2|1), at 23 = 22 + Az the MMP of 
M2 is drawn from /(3|2), and so on. Then, after n steps, 

(S n - So) = (5 c0 Az) 2 v~ 2 ' ( A6 ) 

i=l 

and the distribution of S n is the distribution of the sum of n Inverse Gamma(|,x) variates. Therefore, p n (S n — So = a;) da; 



is 



given by equation (A4), with the obvious notation that 







n 2 x = n 2 (S c0 Az) 2 = S 2 (z n - z ) 2 . (A7) 

Define the formation time tf as the first time that the mass of the MMP becomes less than /Mo, with 0.5 < / < 1. Let 
S{ = <r 2 (/Mo), and let zt denote the redshift associated with the formation time. Suppose that Az is chosen sufficiently small 
that 2f — 20 = nAz, with n an integer. Then, for this MMP generating rule, the formation time distribution is got from 

P(tt < t\MoM) = P(S n -S < Si - So) 

rSf — S 

p n (x) Ax 

= -r-Y M*-*0^ (AS) 
\^/2(Si-S )J 

This is the same as equation (2.23) in Lacey & Cole (1993). 

Suppose instead that the MMP is not chosen according to /(1|0), but it is chosen using iV(l|0) = (Mo/Mi) /(1|0), with 
the requirement that Mo/2 < Mi < Mo. This is the choice we make in the main text. (Recall that, as Az — (zi — 20) — > 0, the 
integral of iV(ljO) -> 1.) Suppose further that at 22 = 21 + Az the MMP of Mi is chosen from JV(2|1), with Mi/2 < M 2 < Mi. 
Then, 

P (S 2 )= f 2 JV(2|l)JV(l|0)d5i, (A9) 

where Smi n = So if S2 < S(Mo/2), and S'min = S(2M2) otherwise. These restrictions on S m i n follow from the requirements 
that, in any step, the mass of the MMP may not decrease to less than half the current mass, and the mass of the MMP can 
never be greater than Mo. Now, if 1S2 < S(Mq/2), then this is the same convolution integral as the split rule by /(1|0) above, 
since the extra factors (Mi /M2) (Mo /Mi) = (M0/M2) may be taken outside the integral. Thus, provided that Ms > Mo/2, 
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p(&) = (j£) /(2|0) = JV(2|0). 



(A10) 



Choose some 0.5 </ < 1. If Mi > M /2 for 1 < i < n, then after n steps, p(S„) = N(S n \So), and 
/•Mo 

P(tf < t\M ,t ) = N(m,Z{\Mo,z )dm, (All) 

where Zf — zo — n Az as before. This expression is the same as equation (2.26) in Lacey & Cole (1993). This demonstrates that, 
for sufficiently small steps in Az, choosing the MMP using 7V(1|0) of equation (^) gives the correct distribution of formation 
times (equation Al). 



A3 Analytic calculation of the distribution of formation masses 

Suppose that the formation time is defined as above for some 0.5 < / < 1. Then the formation mass Alt may be defined in at 
least two ways. The first is to define it as the smallest mass M for which M > /Mo, averaged over all merger histories. This 
is the same definition as that shown in Figure 11 of Lacey & Cole (1993). Alternatively, we could define the formation mass 
as the largest mass M at which M < /Mo, averaged over all histories. Compared to the quantity in the first definition, this 
is the mass of the MMP in the next time step. For either definition, in the appropriate units, the distribution of formation 
masses for a white-noise spectrum (initial slope n = 0) provides an excellent fit to the distribution for other power spectra. 
(Recall that the same is true for the formation time distribution.) This is fortunate, since, when n = 0, then, for our MMP 
generating algorithm, the distribution of formation masses can be calculated analytically. We show this below. 

As above, it is convenient to work in terms of the variable S rather than M. The quantity described by the first definition 
can be computed as follows. The probability that Mf = M is equal to the probability that after step n the mass was 
M n = M > /Mo, and that the mass dropped to M n +i < /Mo on step n + 1, integrated over all allowed values of M„+i, and 
summed over all values of n. The sum over all n allows for all possible formation times. Thus, 

p(M) dM = p{S) dS 

r S(M/2) 

E ' 

allowed n 



AS \ / p(S n +i\S n = S)p(S n = S\S )dS n+1 , (A12) 

I Si 



where Sf — S(fMo). The upper and lower limits to the integral reflect the fact that the allowed values of M„+i are all 
those that are smaller than /Mo but larger than M„/2. Now, the order of the sum and the integral can be rearranged. Since 
p(S n +i\S n = S) is just the one-step split rule, it is independent of n, so it is convenient to do the sum over n first. Since S„ is 
less than St — S(fMo), the previous section showed that, for our algorithm, p(S n \So) dS n = N(S n , z n \So, zo). Moreover, for 
small Az, summing over n is like integrating over all formation times zf. Thus, 



dSj2p(S\S ) ^ J dz t N(S,Zi\S ,zo) 



dS (M /M) IS -So 



(S-So) Az V 2tt5, 



$c0 



(A13) 



This expression is independent of M„ + i. For a white-noise power-spectrum, M/M n+1 = S n +i/S, so the integral over the 
allowed range of S n +i can be done analytically: 

-s(m/2) rjr / Dq - 



J N(S„+i\S„ = S) dS n +i = 2^ — [e- D - —/=-J + C 1 - 2D ) (erf(VDq) - erf (A14) 



where D = {S c0 Az) 2 /2S, and q = S/(S{ - S). Thus, as Az -» 0, so D -> 0, then, provided Dq < 1, 

pf s ) ds -> — dS (2 - — V for n = 0, (A15) 

m ^ y/{s-l)( Si -s) V S J 

where s = S/So, Sf = Si /So, and m = M/Mo- In this limit, the dependence on Az has dropped out. 

The quantity described by the second definition can be computed similarly. The probability that Mf = M is equal to the 
probability that after step n the mass was M n > /Mo, and that the mass dropped to M < /Mo on step n + 1, integrated 
over all allowed values of M n , and summed over all choices of n. That is, 

p{M) dM = p(S) dS 

= dS / p(S n+ i = S\S n )p(S n )dS n , (A16) 

allowed „JSmin 
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Figure A2. Di stribu tion of formation masses using our MMP generating algorithm. The thick solid curve that is lowest on the left 
shows equation (A19) with / = 0.5 and M = fMo- Curves for other values of / are similar so we have not shown them. The other thick 
curves show equation ( |A15| ) with / = 0.5, 0.7 and 0.9, and M = Mq. The thin curves show the corresponding quantities computed using 
our Monte-Carlo algorithm for a range of values for the spectral slope n (0, —1 and —2), the final mass Mo (0.1, 1, and 10Af„), and the 
time step Az (0.01—0.1). All the curves are very similar, and they are well fit by the analytic (white-noise) formulae. 



where S'min = S(2M) if 2M < Mo, and S m in = So otherwise because, for our algorithm, Mf — M must lie in the range 
/M /2 < M < /Mo. This means that M n lies in the range M < M„ < 2M, unless 2M > M , in which case the maximum 
allowed value of M n is simply Mo- The sum is over all formation times, so it is over all n > 0. Now, as before, p(S n ) dS„ = 
N(S„, z n \So, zo), and p(S„+i = S\S n ) is just the one step split rule, so it is independent of n. If we rearrange the order of the 
sum and the integral above, then, for small Az, summing over n is like integrating over all formation times zi, so the integral 
is the same as before: 



dS„ ^2p(S„) 



dS n 
Az 



dz{ N{S n ,Z{\So,z ) 



dS n (Mo/Mn) S n -So 



2<, 



(5„ - So) Az 

Using this expression and then doing the integral over the allowed range of S n implies that 



p(S)dS - 



dS 



/Mo 



(A17) 



(A18) 



where 



A 



(Sep Az) 2 
2(5 - So) ' 



B 



[ Si-S \ 
\S-Sf ) 



and C = max ( 0, 1 



For small Az, the error function terms become 2^/[A/tt)(VB — VC) provided (S — Si) /Si is not small. In this approximation, 
the formation mass distribution is 



p(s) ds 



(Mi) 1 

\ M ) 7r(s - So) V s 



SO 



ds 



(Ml) 1 

V M J tt(s - So) 



so 



s - 1 



so 

■'"'mm 



if 2M > M , 



ds if M > 2M > /Mo, 



(A19) 
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where Si = Si/Sf. Notice that, in this limit, the dependence on Az has dropped out. 

As with the distribution of formation times, this is a well-behaved probability distribution only for white-noise, where 
Si = Sj/iSf = f Mo/Mi. However, as with the formation time distribution, this distribution depends only very weakly on power 
spectrum, so the white-noise expression provides a good fit to the distribution for other power-spectra. It also turns out that 
this distribution depends only weakly on /, so, in Figure A2, only the curve for / = 0.5 is shown. 

Figure A2 shows both formation mass distributions for a range of initial power-spectra, final Mos, and choices for /. The 
a>axis shows Mf/Mo for equation (A15), and Alt / (f Mo) for equation (A18). The thick solid curves show the analytic formulae 
with / = 0.5, 0.7 and 0.9, for M/M* = 0.1, 1, and 10, and n = 0. The thinner lines show the associated distributions obtained 
using our MMP generating code for n = (solid), n = — 1 (dashed), and n = —2 (dotted). The curves for the different 
power-spectra are virtually indistinguishable, and they are all well fit by our analytical formulae. This, and the fact that our 
code also generates the correct distribution of formation times (Fig. 1), shows that our MMP generating code works. 
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